rm(list=ls(all=T))
options(stringsasFactors = FALSE)

p_needed <- c("dplyr", "stargazer", "reshape2", "brms","ggplot2")
lapply(p_needed, require, character.only = TRUE)

dat <- readRDS("cleandata_replication.Rds")

## Figure 3: Predicted probabilities labor market risks ----
econ <- readRDS("fit/econ30.Rds")
posterior <- posterior_samples(econ)
posterior <- posterior[,grep("b_", colnames(posterior))]

econfix <- matrix(nrow=100, ncol=85, rep(0, 85*100))
colnames(econfix) <- row.names(fixef(econ))[grep("factorcyear", row.names(fixef(econ)))]
grep("SouthAfrica2008", colnames(econfix))
econfix[,67] <- 1

#adjust for EA?
#person 1, no formal schooling
e.ct1 <- cbind(1, #intercept
               0, #educ = No formal schooling
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100), #night lights
               37, #age 
               0, #male
               1, #urban
               0, #unemployed
               0, #ethnic grievances
               6, #media consumption
               0, #conflict
               mean(dat$ea_service, na.rm = T),
               mean(dat$ea_infrastructure, na.rm = T),
               econfix,
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100)*0 #interaction
)
#person 2,  secondary education
e.ct2 <- cbind(1, #intercept
               5, #educ = secondary schooling
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100), #night lights
               37, #age 
               0, #male
               1, #urban
               0, #unemployed
               0, #ethnic grievances
               6, #media consumption
               0, #conflict
               mean(dat$ea_service, na.rm = T),
               mean(dat$ea_infrastructure, na.rm = T),
               econfix,
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100)*5 #interaction
)
#person 3, university degree
e.ct3 <- cbind(1, #intercept
               8, #educ = post-grad
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100), #night lights
               37, #age 
               0, #male
               1, #urban
               0, #unemployed
               0, #ethnic grievances
               6, #media consumption
               0, #conflict
               mean(dat$ea_service, na.rm = T),
               mean(dat$ea_infrastructure, na.rm = T),
               econfix,
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100)*8 #interaction
)
posterior <- as.matrix(posterior)

y.ct1 <- posterior %*% t(e.ct1)
y.ct2 <- posterior %*% t(e.ct2)
y.ct3 <- posterior %*% t(e.ct3)

y.ct1 <- exp(y.ct1)/(1 + exp(y.ct1))
y.ct2 <- exp(y.ct2)/(1 + exp(y.ct2))
y.ct3 <- exp(y.ct3)/(1 + exp(y.ct3))

#sphaghetti plot
cib <- 0.05
lb <- round((dim(posterior)[1] * cib)/2)
ub <- dim(posterior)[1] - lb

rr <- sample(lb:ub, 300)
if(!lb%in%rr) rr <- c(rr,lb)
if(!ub%in%rr) rr <- c(rr,ub)

dat1 <- cbind(Light=seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), 
                       length.out = 100), 
             colMeans(y.ct1), colMeans(y.ct2), colMeans(y.ct3)) %>% 
  as.data.frame 

dat1 <- melt(dat1, id.vars="Light")
dat1$variable <- as.character(dat1$variable)
dat1$variable[dat1$variable=="Light"] <- "Light"
dat1$variable[dat1$variable=="V2"] <- "No Formal Schooling"
dat1$variable[dat1$variable=="V3"] <- "Secondary School"
dat1$variable[dat1$variable=="V4"] <- "University Degree"

dat1$ind <- rep(1:100, 3)

e.int <- ggplot(dat1, aes(x=Light, y=value, color=variable)) +
  facet_wrap(~variable)+
  scale_color_manual(values=c("black", "black", "black"))+
  geom_line() +
  theme_bw() +
  labs(y="Pr(Labor market risks)", x="Local Development (Logged Night Light)")+
  theme(legend.position = "none") +
  scale_y_continuous(limits = c(0,1))

panel.grid.major = element_blank()

options(expressions=10000)
for(i in rr){ #lb:ub
  df <- cbind(Light=seq(min(dat$light_30_log, na.rm=T), 
                        max(dat$light_30_log, na.rm=T), 
                        length.out = 100), 
              y.ct1[i,], 
              y.ct2[i,], y.ct3[i,]) %>% as.data.frame
  dat2 <- melt(df, id.vars="Light")
  dat2$variable <- as.character(dat2$variable)
  dat2$variable[dat2$variable=="Light"] <- "Light"
  dat2$variable[dat2$variable=="V2"] <- "No Formal Schooling"
  dat2$variable[dat2$variable=="V3"] <- "Secondary School"
  dat2$variable[dat2$variable=="V4"] <- "University Degree"
  dat2$ind <- rep(1:100, 3)
  
  e.int <- e.int + 
    geom_line(data=dat2, aes(x=Light, y=value, color=variable), size=1, alpha=0.02)
  print(i) 
}

e.int


## Figure 4: Predicted probabilities living conditions ----
liv <- readRDS("fit/liv30.Rds")
posterior <- posterior_samples(liv)
posterior <- posterior[,grep("b_", colnames(posterior))]
intercepts <- posterior[,grep("b_Intercept", colnames(posterior))]
posterior <- posterior[,-c(1:4)]

livfix <- matrix(nrow=100, ncol=85, rep(0, 85*100))
colnames(livfix) <- row.names(fixef(liv))[grep("factorcyear", row.names(fixef(liv)))]
grep("SouthAfrica2008", colnames(livfix))
livfix[,67] <- 1

#adjust for EA?
#person 1, no formal schooling
e.ct1 <- cbind(0, #educ = No formal schooling
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100), #night lights
               37, #age 
               0, #male
               1, #urban
               0, #unemployed
               0, #ethnic grievances
               6, #media consumption
               0, #conflict
               mean(dat$ea_service, na.rm = T),
               mean(dat$ea_infrastructure, na.rm = T),
               livfix,
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100)*0 #interaction
)
#person 2,  secondary education
e.ct2 <- cbind(5, #educ = secondary schooling
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100), #night lights
               37, #age 
               0, #male
               1, #urban
               0, #unemployed
               0, #ethnic grievances
               6, #media consumption
               0, #conflict
               mean(dat$ea_service, na.rm = T),
               mean(dat$ea_infrastructure, na.rm = T),
               livfix,
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100)*5 #interaction
)
#person 3, university degree
e.ct3 <- cbind(8, #educ = post-grad
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100), #night lights
               37, #age 
               0, #male
               1, #urban
               0, #unemployed
               0, #ethnic grievances
               6, #media consumption
               0, #conflict
               mean(dat$ea_service, na.rm = T),
               mean(dat$ea_infrastructure, na.rm = T),
               livfix,
               seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), length.out = 100)*8 #interaction
)
posterior <- as.matrix(posterior)

y.ct1 <- intercepts[,3] - posterior %*% t(e.ct1)
y.ct2 <- intercepts[,3] - posterior %*% t(e.ct2)
y.ct3 <- intercepts[,3] - posterior %*% t(e.ct3)

y.ct1 <- 1/(1 + exp(y.ct1))
y.ct2 <- 1/(1 + exp(y.ct2))
y.ct3 <- 1/(1 + exp(y.ct3))

cib <- 0.05
lb <- round((dim(posterior)[1] * cib)/2)
ub <- dim(posterior)[1] - lb

rr <- sample(lb:ub, 300)
if(!lb%in%rr) rr <- c(rr,lb)
if(!ub%in%rr) rr <- c(rr,ub)

dat1 <- cbind(Light=seq(min(dat$light_30_log, na.rm=T), max(dat$light_30_log, na.rm=T), 
                       length.out = 100), 
             colMeans(y.ct1), colMeans(y.ct2), colMeans(y.ct3)) %>% 
  as.data.frame 

dat1 <- melt(dat1, id.vars="Light")
dat1$variable <- as.character(dat1$variable)
dat1$variable[dat1$variable=="Light"] <- "Light"
dat1$variable[dat1$variable=="V2"] <- "No Formal Schooling"
dat1$variable[dat1$variable=="V3"] <- "Secondary School"
dat1$variable[dat1$variable=="V4"] <- "University Degree"

dat1$ind <- rep(1:100, 3)

e.int <- ggplot(dat1, aes(x=Light, y=value, color=variable)) +
  facet_wrap(~variable)+
  scale_color_manual(values=c("black", "black", "black"))+
  geom_line() +
  theme_bw() +
  labs(y="", x="Local Development (Logged Night Light)")+
  theme(legend.position = "none") +
  scale_y_continuous(limits = c(0,1))

options(expressions=10000)
for(i in rr){ #lb:ub
  df <- cbind(Light=seq(min(dat$light_30_log, na.rm=T), 
                        max(dat$light_30_log, na.rm=T), 
                        length.out = 100), 
              y.ct1[i,], 
              y.ct2[i,], y.ct3[i,]) %>% as.data.frame
  dat2 <- melt(df, id.vars="Light")
  dat2$variable <- as.character(dat2$variable)
  dat2$variable[dat2$variable=="Light"] <- "Light"
  dat2$variable[dat2$variable=="V2"] <- "No Formal Schooling"
  dat2$variable[dat2$variable=="V3"] <- "Secondary School"
  dat2$variable[dat2$variable=="V4"] <- "University Degree"
  dat2$ind <- rep(1:100, 3)
  
  e.int <- e.int + 
    geom_line(data=dat2, aes(x=Light, y=value, color=variable), size=1, alpha=0.02)
  print(i) 
}
e.int

## Figure 5: Predicted probabilities labor market risks (growth) ----
econdyn <- readRDS("fit/econ3ch.Rds")
posterior <- posterior_samples(econdyn)
posterior <- posterior[,grep("b_", colnames(posterior))]

fe <- matrix(nrow=100, ncol=78, rep(0, 78*100))
colnames(fe) <- row.names(fixef(econdyn))[grep("factorcyear", row.names(fixef(econdyn)))]

#find year and country
grep("SouthAfrica2008", colnames(fe))
fe[,61] <- 1

#adjust for EA?
#person 1, no formal schooling
e.ct1 <- cbind(1, #intercept
               0, #educ = No formal schooling
               seq(min(dat$change3_30_log, na.rm=T), max(dat$change3_30_log, na.rm=T), length.out = 100), #night lights
               37, #age 
               0, #male
               1, #urban
               0, #unemployed
               0, #ethnic grievances
               6, #media consumption
               0, #conflict
               mean(dat$ea_service, na.rm = T),
               mean(dat$ea_infrastructure, na.rm = T),
               fe,
               seq(min(dat$change3_30_log, na.rm=T), max(dat$change3_30_log, na.rm=T), length.out = 100)*0 #interaction
)
#person 2,  secondary education
e.ct2 <- cbind(1, #intercept
               5, #educ = secondary schooling
               seq(min(dat$change3_30_log, na.rm=T), max(dat$change3_30_log, na.rm=T), length.out = 100), #night lights
               37, #age 
               0, #male
               1, #urban
               0, #unemployed
               0, #ethnic grievances
               6, #media consumption
               0, #conflict
               mean(dat$ea_service, na.rm = T),
               mean(dat$ea_infrastructure, na.rm = T),
               fe,
               seq(min(dat$change3_30_log, na.rm=T), max(dat$change3_30_log, na.rm=T), length.out = 100)*5 #interaction
)
#person 3, university degree
e.ct3 <- cbind(1, #intercept
               8, #educ = post-grad
               seq(min(dat$change3_30_log, na.rm=T), max(dat$change3_30_log, na.rm=T), length.out = 100), #night lights
               37, #age 
               0, #male
               1, #urban
               0, #unemployed
               0, #ethnic grievances
               6, #media consumption
               0, #conflict
               mean(dat$ea_service, na.rm = T),
               mean(dat$ea_infrastructure, na.rm = T),
               fe,
               seq(min(dat$change3_30_log, na.rm=T), max(dat$change3_30_log, na.rm=T), length.out = 100)*8 #interaction
)
posterior <- as.matrix(posterior)

y.ct1 <- posterior %*% t(e.ct1)
y.ct2 <- posterior %*% t(e.ct2)
y.ct3 <- posterior %*% t(e.ct3)

y.ct1 <- exp(y.ct1)/(1 + exp(y.ct1))
y.ct2 <- exp(y.ct2)/(1 + exp(y.ct2))
y.ct3 <- exp(y.ct3)/(1 + exp(y.ct3))

cib <- 0.05
lb <- round((dim(posterior)[1] * cib)/2)
ub <- dim(posterior)[1] - lb

rr <- sample(lb:ub, 300)
if(!lb%in%rr) rr <- c(rr,lb)
if(!ub%in%rr) rr <- c(rr,ub)

dat1 <- cbind(Light=seq(min(dat$change3_30_log, na.rm=T), max(dat$change3_30_log, na.rm=T), 
                       length.out = 100), 
             colMeans(y.ct1), colMeans(y.ct2), colMeans(y.ct3)) %>% 
  as.data.frame 

dat1 <- melt(dat1, id.vars="Light")
dat1$variable <- as.character(dat1$variable)
dat1$variable[dat1$variable=="Light"] <- "Light"
dat1$variable[dat1$variable=="V2"] <- "No Formal Schooling"
dat1$variable[dat1$variable=="V3"] <- "Secondary School"
dat1$variable[dat1$variable=="V4"] <- "University Degree"

dat1$ind <- rep(1:100, 3)

e.int <- ggplot(dat1, aes(x=Light, y=value, color=variable)) +
  facet_wrap(~variable)+
  scale_color_manual(values=c("black", "black", "black"))+
  geom_line() +
  theme_bw() +
  labs(y="", x="Change in Local Development (Logged Night Light)")+
  theme(legend.position = "none") +
  scale_y_continuous(limits = c(0,1))

options(expressions=10000)
for(i in rr){ #lb:ub
  df <- cbind(Light=seq(min(dat$change3_30_log, na.rm=T), 
                        max(dat$change3_30_log, na.rm=T), 
                        length.out = 100), 
              y.ct1[i,], 
              y.ct2[i,], y.ct3[i,]) %>% as.data.frame
  dat2 <- melt(df, id.vars="Light")
  dat2$variable <- as.character(dat2$variable)
  dat2$variable[dat2$variable=="Light"] <- "Light"
  dat2$variable[dat2$variable=="V2"] <- "No Formal Schooling"
  dat2$variable[dat2$variable=="V3"] <- "Secondary School"
  dat2$variable[dat2$variable=="V4"] <- "University Degree"
  dat2$ind <- rep(1:100, 3)
  
  e.int <- e.int + 
    geom_line(data=dat2, aes(x=Light, y=value, color=variable), size=1, alpha=0.02)
  print(i) 
}
e.int


# Posterior distribution main models----
# Figure 6: Living conditions
liv <- readRDS("fit/liv30.Rds")
test <- as.matrix(coda::as.mcmc(liv))
test <- as.data.frame(test)
test <- test %>%
  dplyr::rename(Education=b_educ, Age = b_age,Urban = b_urban,Unemployment = b_unempl, Conflict = b_conflict)
names(test)[names(test) == "b_light_30_log"] <- "Local Development" 
names(test)[names(test) == "b_ethfair"] <- "Ethnic Grievances"
names(test)[names(test) == "b_educ:light_30_log"] <- "Education * Local Development"
names(test)[names(test) == "b_ea_service"] <- "Services (EA)"
names(test)[names(test) == "b_ea_infrastructure"] <- "Infrastructure (EA)"
names(test)[names(test) == "b_media"] <- "Media Consumption"

bayesplot::color_scheme_set("gray")
bayesplot::mcmc_areas(test, pars = c("Education","Local Development","Education * Local Development","Age","Urban","Unemployment","Media Consumption",
                                     "Ethnic Grievances","Services (EA)","Infrastructure (EA)", "Conflict"),
                      prob = 0.5) +
  labs(title = "", subtitle = "") +
  geom_vline(xintercept = 0, size = 0.5, linetype = "dashed") +
  theme_bw()

# Figure 7: Labor market risks
econ <- readRDS("fit/econ30.Rds")
test <- as.matrix(coda::as.mcmc(econ))
test <- as.data.frame(test)
test <- test %>%
  dplyr::rename(Education=b_educ, Age = b_age,Urban = b_urban,Unemployment = b_unempl, Conflict = b_conflict)
names(test)[names(test) == "b_light_30_log"] <- "Local Development" 
names(test)[names(test) == "b_ethfair"] <- "Ethnic Grievances"
names(test)[names(test) == "b_educ:light_30_log"] <- "Education * Local Development"
names(test)[names(test) == "b_ea_service"] <- "Services (EA)"
names(test)[names(test) == "b_ea_infrastructure"] <- "Infrastructure (EA)"
names(test)[names(test) == "b_media"] <- "Media Consumption"

# Coef Distribution Plots
bayesplot::color_scheme_set("gray")
bayesplot::mcmc_areas(test, pars = c("Education","Local Development","Education * Local Development","Age","Urban","Unemployment","Ethnic Grievances",
                                     "Media Consumption","Services (EA)","Infrastructure (EA)", "Conflict"),
                      prob = 0.5) +
  labs(title = "", subtitle = "") +
  geom_vline(xintercept = 0, size = 0.5, linetype = "dashed") +
  theme_bw()

